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We propose a scheme for the initiahzation of a quantum computer based on neutral atoms trapped 
in an optical lattice with large lattice constant. Our focus is the development of a compacting 
scheme to prepare a perfect optical lattice of simple orthorhombic structure with unit occupancy. 
Compacting is accomplished by sequential application of two types of operations: a flip operator 
that changes the internal state of the atoms, and a shift operator that moves them along the lattice 
principal axis. We propose physical mechanisms for realization of these operations and we study 
the effects of motional heating of the atoms. We carry out an analysis of the complexity of the 
compacting scheme and show that it scales linearly with the number of lattice sites per row of the 
lattice, thus showing good scaling behavior with the size of the quantum computer. 

I. INTRODUCTION 

The present work is motivated by the possibility of using neutral (alkali) atoms trapped in an optical lattice for 
quantum computing. Q, |^ |^ We focus on a scalable quantum computing architecture where the number of qubits 
can in principle be increased using only polynomially more physical resources. We consider one-dimensional (1-D), 
orthorhombic (in the near cubic arrangement) two-dimensional (2-D) and three-dimensional (3-D) lattices, although 
generalization of the proposed scheme to other structures is possible. The lattice is characterized by a large lattice 
constant a which allows each lattice qubit to be individually addressed during the computation. The total number of 
sites N is defined as the number of sites within the minimal compact volume enclosing the region of the lattice that 
is populated by atoms. Assuming that the lattice sites are either occupied by a single atom or are vacant, we define a 
filling factor / — Nocc/N, where Nqcc is the number of (singly) occupied sites. Under the assumption that all lattice 
sites are either occupied with a single atom or empty (i.e. with a miiform distribution over the lattice), the filling 
factor is the same as the site occupancy probability Pocc- 

We investigate here the preparation of the system for quantum computation, i.e. its initialization as a quantum 
computer. Our key objective is a perfect optical lattice where each lattice site is occupied with a single atom in 
a specific internal state and the motional ground state. We present a feasible and scalable scheme for compacting 
the optical lattice to achieve this state with 100% fidelity. Our scheme removes the vacant sites to the edge of the 
initial lattice {Ni, fi) and thus creates a smaller but perfect final lattice {Nf,ff), where Nf < Ni and ff — 1 > fi. 
This compacting procedure can be used for initialization, or for reinitialization of a quantum computer before a 
new computation. Other applications may also be found in the construction of fault-tolerant procedures for the 
computation itself. 

A motivation for developing the compacting procedure is the conjecture that randomly occuring imperfections, 
even if their locations are known, cause bottlenecks in quantum information flow which can not be avoided when 
the computer size, i.e. the number of lattice qubits, scales up. In fact, we maintain that the probability of finding 
a "good" sublattice is exponentially small for any constant filling factor. This is because the probability that there 
will be 'insurmountable' blocks of defects (gaps) in any chosen sublattice increases rapidly with its size. The site 
percolation threshold Pc reaches the following values: 1 (1-D), 0.59 (2-D) and 0.31 (3-D) Thus the initial filling 

factor fi = 0.5 exceeds the percolation threshold only in a three dimensional lattice. Even in this case, the condition 
fi > Pc only implies a non-zero probability that distant qubits may be connected. It does not ensure that they are 
connected via independent routes, nor even that they are connected at all. In addition, the mapping of a quantum 
algorithm onto an imperfect lattice may be a hard classical computational problem (even NP hard). Based on these 
considerations, an imperfect lattice structure seems to pose a stumbling block in the realization of a scalable quantum 
computer. 

We briefly list other possible approaches to preparing an optical lattice with single occupancy at each site. The 
quantum phase transition from the supcrfluid to Mott insulator phase ^ 0| , recently observed in (sj , can prepare a 
singly occupied optical lattice when both the lattice constant and the well depth are rather small (a w 0.5 fim, Uq ~ 1 
/iK). The quantum critical point is given by the ratio of on-site atomic interaction energy U and the tunneling strength 
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T as U /T = 10.6 d, where d is the dimensionahty of the lattice With increasing values of lattice constant or 

well depth, the tunneling strength diminishes exponentially, bringing the system deeply into the Mott insulator phase. 
Here, the band structure of the lattice potential energy spectrum approaches its discrete limit, and is suitable for 
robust quantum computation with independent qubits. This approach can prepare a singly occupied optical lattice 
of small lattice constant with more than 90% of fidelity. However, it is difficult to apply it to an addressable optical 
lattice, where the large lattice constant sets an unrealistically long timescale for adiabatically changing the lattice. In 
addition, even if the phase transition could be established, the final singly occupied optical lattice would be unlikely 
to reach the desired level of perfection, and subsequent compacting would be required. Another recent proposal for 
preparing a nearly perfect optical lattice is based on use of the dipole interaction between atoms excited to Rydberg 
states Here, quantum state manipulations exploiting the dipole blockade effect can lead to a lattice with a 

significantly fewer vacant sites |lOj |. Neither of the two previously presented approaches can currently guarantee a 
perfect lattice where each site is occupied by exactly one atom. An alternative proposal for preparing a perfect optical 
lattice has recently been proposed. It involves adiabatic loading of one optical lattice from another that has one or 
more atoms at every site [llj . 



II. PHYSICAL SYSTEM 



An optical lattice with individually addressable sites is a particularly promising candidate for implementation of 
quantum computation. This important requirement is satisfied by having a large lattice constant (a 5 //to). A 1-D 
optical lattice can easily be realized by interfering two linearly polarized laser beams |12| . Higher dimensional lattices 
can be implemented using two (2-D) or three (3-D) pairs of laser beams in a mutually perpendicular arrangement. 
Undesired interference between sublattices can be eliminated by using slightly different frequencies for each of the 
beam pairs. The resulting lattice possesses a simple orthorhombic structure. Such an optical lattice with a « 5 /xm 
can be made with C02-laser beams or with blue-detuned light. The blue detuned standing waves consist of two 
beams propagating at a shallow angle Of, with respect to each other, giving a lattice spacing a = \/2sin{9ii/2). 

A pair of counterpropagating (along the z axis) linearly polarized laser beams of identical wavelength A gives rise 
to a superposition of standing waves of opposite circular polarization e± . The electric field is given by: 



E{z) = V2Ea[e'^/'^cos{kz - 0/2)e^ - e'^^'^cosikz + 9/2)e+]. (1) 

Here 6 is the relative angle between the linear polarization vectors of both beams, k = 27r/A, and Et^ is the single beam 
field amplitude. In the absence of any additional external field, the resulting 1-D periodic lattice potential depends 
on the magnetic hyperfine sublevel [l|. It can be characterized by the following relation 

U{z) = ^cos{9)cos{2kz) + ^sin{e)sin{2kz), (2) 

where C/o — jOtE"^ and Ui — ^aE^^-^ describe the well depths of the potential at 6* = and 6 ^ tt respectively, 
a is the characteristic polarizability of a given transition, F is the total angular momentum of the relevant atomic 
hyperfine level and m p is the magnetic hyperfine sublevel. Note that the potentials for all magnetic hyperfine sublevels 
coincide for = 0. 

One-dimensional as well as multi-dimensional orthorhombic arrangements with approximately 20 sites per dimension 
and 150 /iK depth are readily achievable. For 1-D, 2-D and 3-D lattice potentials and a filling factor of 0.5, this 
results in 10, 200, and 4000 atoms, respectively. We focus on an optical lattice filled with atoms of ^^'^Cs, although 
the proposed scheme is applicable to other alkali atom. The 65^6*1/2 electronic ground state of ^^'^Cs consists of two 
hyperfine levels of total angular momentum F = 3 and F = 4, with energy splitting AF = Ef=4 — Ep=^ = 9.1926 
GHz. 



III. LOADING AND COOLING OF OPTICAL LATTICE 

A magneto-optical trap |13| can be used to load well-spaced far-off-resonant optical lattice sites with man y at oms 
(N=10-100) Atoms are lost during laser cooling in the lattice, in pairs, due to photon assisted collisions [islllq . 
Eventually, the sites initially occupied by an odd number of atoms become singly occupied, and initially even occupied 
sites are left empty. The filling factor, /i, that results from this combined loading and cooling process approaches 
a maximum of 0.5. This sort of loading of a far-off-resonant optical lattice with Cs atoms, to 44% filling, has been 
demonstrated experimentally with a small lattice constant (~ 0.5 /im). 
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FIG. 1: Compression scheme, a) All atoms of the lattice are in a single internal state and are shifted by +a/2 as a result of 
rotation of the relative polarization angle 9 from to vr. b) The state of the mobile atoms is then flipped, e.g. from rap = +1 
(filled circles) to mp = — 1. This is accompanied by a change in shift direction, denoted by the arrows, c) Rotation of the angle 
6 back from tt to moves the = +1 atom back to its original position and the rap = — 1 atom forward to the next lattice 
site, providing the desired compacting of the lattice structure, d) The mobile atoms are then flipped back to the original state. 
The net effect is a shift of the right hand atom by one site to the left, or equivalently, of the corresponding vacancy by one site 
to the right. 



Trapped atoms can then be brought to the vibrational ground state of their lattice sites using Raman sideband 
coolling 0, 0, IT^ I20I l2l|. This procedure has been experimentally demonstrated in ID [HiiJ, 2D and 3D 
[13, 1201 optical lattices. In [23, the 3D ground state was populated by up to 55% of atoms, corresponding to over 80% 
in each dimension. It was shown that the cooling was limited by the rescattering of cooling photons. This mechanism 
is dramatically reduced at the low densities associated with large lattice constants. Even with a 5 /xm lattice constant, 
a 150 /zK depth lattice leaves the atoms well within the Lamb-Dicke limit required for efficient sideband cooling. It 
is reasonable to expect nearly 100% vibrational ground state population after Raman sideband cooling. 

With a 150 deep, 5 /xm spaced lattice, the vibrational level spacing is well below the expected minimum 
polarization gradient cooling temperature of 2 /iK, so the polarization gradient cooling limit should be obtained 
[l^. With a temperature so small compared to the lattice depth, there will be negligible site hopping during cooling. 
Therefore the atoms can be imaged during the cooling for as long as needed. ID CO2 lattices have recently been 
imaged |12| . In a similar way, it should be possible to use ~0.6 numerical aperture optics to image successive 2D 
planes, and thus determine the site occupancy of a 3D lattice. Since the localized atoms in the 3D lattice will scatter 
light coherently, there will be significant interference between light scattered from atoms in the image plane and light 
from non-imaged planes. Information can be obtained from such signals, easier so if the lattice spacing can be made 
an integer multiple of image wavelengths, as is possible in a blue-detuned lattice. The total intensity of light from 
each non-imaged plane is only about 2% of the intensity from the image plane, making the signal to background from 
a 4000 atom 3D lattice ~3. When the lattice occupation is perfectly known, it can then be compacted into a perfectly 
occupied lattice. 



IV. COMPACTING PROCEDURE 



We now present a scheme for converting an imperfect optical lattice, where each site is either empty or occupied 
by a single atom, into a smaller lattice where each site is occupied by a single atom. Our scheme exploits the ability 
to move a subset of trapped atoms to fill vacant sites, thereby removing vacancies to the edge of the lattice potential. 
Compacting a lattice is equivalent to sequentially removing vacant sites to lattice edges. 

We define two elementary operations that are sufficient to compact the lattice. They are (1) the flip operator, F/, 
which toggles the state of an atom at position / = (i, j, k) between two levels, and (2) the shift operator, Sij, which 
moves an atom from position / to a neighboring position J = k'), where only one of the coordinates of J differs 
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from / by one. 

Compacting results from translating some of the atoms, which we call mobile atoms, to fill vacant sites, while 
the rest, which we call stationary atoms, remain fixed within their lattice structure. We refer to the state of a 
mobile atom as the mobile state, and that of a stationary atom as the storage state. These may for instance be 
\Cs 6Si/2,F = 3,mp = —1) and \Cs 6Si/2,F = 3,mp — +1) respectively. The mobile atoms are initially selected 
based on the map of lattice occupancy obtained from the imaging process. The compacting process starts with all 
atoms in the storage state. Atoms are moved by rotating the linear polarization of one of the lattice beams by angle 
9 relative to the counterpropagating beam, according to the following steps: 

1. for ^ = — > TT, all atoms are initially in the storage state and all move by —a/2, where a is the lattice constant 
and the sign indicates direction; formally, this operation keeps the lattice occupation structure invariant and 
hence acts on it as identity; 

2. at 6* = TT, the mobile atoms are selectively flipped into the mobile state, while the states of stationary atoms 
remain unchanged; 

3. for 9 = TT 0, the stationary atoms are shifted back to their original positions (+a/2) while the mobile atoms 
are shifted forward to the next vacant lattice site (—a/2); 

4. the mobile atoms are selectively flipped back to the storage state. 

After step 3, the lattice occupation structure is now different. All atoms that were selected to be mobile have 
been shifted to vacant sites. Vacancies have been displaced to the edge of the lattice. The procedure is repeated 
until all vacant sites have been moved to edge of the lattice. This scheme is summarized in Fig. ^ For practical 
implementation, the scheme can be simplified by not requiring that storage atoms end each operation in the same 
position. As the lattice is compacted to its center, translations of the storage atoms will tend cancel each other out. 

In the next two sections we consider the physics of each elementary step in detail in order to place bounds on the 
extent of undesired heating. 

V. FLIP OPERATION 
A. Achieving site selectivity 

The compacting scheme requires that we can make an atom at a single site undergo a flip transition, while none 
of its neighbors make the transition. We propose two distinct ways to accomplish this feat. The first, and more 
novel, approach uses an independent "addressing" beam, which is a far-off-resonant, circularly polarized laser beam 
tightly focused on the atom to be flipped. The circularly polarized beam shifts the stationary and mobile atomic 
states differently, so that the resonance transition between them is shifted. The addressing beam will have non-zero 
intensity at non-target atoms, especially those than lie along the addressing beam axis. But as long as the Rayleigh 
range is reasonably much smaller than the lattice spacing, the resonance frequency shift will be much smaller at 
non-target sites. In the presence of the addressing beam, a pulse from a spatially homogeneous source can be used to 
flip the target atom. The requirement is that the pulse time be long compared to the inverse of the difference between 
the resonance frequency shifts of the target and the nearest non-target atoms. In that case, non-target atoms will 
be far enough off resonance that they will not be flipped. These flip pulses can be driven either by direct microwave 
excitation or by co-propagating stimulated Raman beams. The required parameters of an addressing beam are easily 
attained. For instance, a 2 /itW laser beam at 877 nm, focused to a waist of 1.2 fim, gives a 1 MHz relative Stark 
shift, 2.5 times larger than that of the nearest neighbor. 

Another way to get site selectivity is to drive an off-resonant stimulated Raman transition using tightly focused 
perpendicular laser beams, so that only atoms at the target site see appreciable intensity from both beams. An 
advantage to this approach is that, since there is no need to frequency resolve the transitions at different sites, the 
flip can be accomplished much more quickly. We will discuss some disadvantages in the following subsections. 

B. Driving flip transitions 

All of the various ways to make site-selective flips have the option of driving a 7r-pulse or using adiabatic fast 
passage. With 7r-pulses, the best approach is probably to tailor the pulse shape, as with a Blackman pulse, in order 
to minimize off-resonant excitation. The option also exists to use a square pulse, but it requires to make sure that all 
non-target atoms lie close to the first minimum of the resulting sine function |22| | . Either way, 7r-pulses require very 
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FIG. 2: The three-level quantum system coupled by two optical fields in a A-configuration. flij is the Rabi frequency for 
coupling of the states i and j, Aij is the detuning from a resonance cuij — Ei — Ej /h. 



stable and repeatable field intensity. In the case where Raman beams are used without an addressing beam, stable 
and repeatable centering of the tightly focused beams on the target atom would be a significant challenge for making 
a reliable flip. 

The flip operation can also be implemented using rapid adiabatic passage. Consider a three level optical system in 
a A-configuration (see Fig. ^ where two ground states |0) and |2) are coupled indirectly via an intermediate state 
using a pair of laser fields. Adiabatic passage occurs when the initial and target levels coupled off-resonantly through 
the intermediate state sweep through their dressed-state resonance as result of interaction with a line arly frequency 
chirped field (for a description of the chirped field and its interaction with matter see for instance |23.l24|'). 

The atom-field Hamiltonian is (in the interaction representation) 
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H(t) = Ho -f Hi(t) = > hio,\i){i\ + [VFoi(i)(7+ + h.c] + [VFi2(t)(7+ + h.c] (3) 



1=0 



where crfj = \j){i\ and Wij{t) is the time-dependent coupling strength given by the external optical field. Canonical 
transformation of this system with coupling fields detuned from resonance with the intermediate state results in an 
effective two-level system [2^ : 

Heff(t) — Ho + Hstark + Hi^eff (i) 

2 

i=Q 

\Wo^it)\' , , \W,,it)\^ , 

+ T f^01 + T <^12 

^01 ^12 

+ Wo2a+ -f WoVo2 (4) 

where W02 = '^oiWi^ (— L- — — L.) jg the effective coupling strength between the states |0) and |2), and Ay is the 
detuning of the field from the atomic transition between levels i and j . Application of a linearly chirped field for one 
of the transitions ^ 1 or 2 — > 1 results in a robust and complete flip operation • Such an adiabatic rapid passage 
can also be done with microwave radiation. 
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C. Heating due to flip transitions 

Site selective flip operations can cause heating in several distinct ways. The first is due to the impulse that both 
target and non-target atoms can feel during the process of turning on the spatially inhomogeneous beams, including 
either the addressing beam or the tightly focused stimulated Raman beams. However, when the frequency of these 
beams is tuned between the first excited state fine structure levels {QP1/2 and 6P3/2 in Cs) the ac Stark shifts due to 
the two levels are opposite. For any individual ground state sublevel, there is a magic frequency where the two Stark 
shifts cancel. To avoid this heating, one simply needs to use the magic frequency for the storage hyperfine sublevel, 
for instance, 877 nm for the F = 3, mp = 1 hyperfine sublevel. The difference in frequency between the Raman 
beams is negligible on the scale needed to avoid this heating effect. 

A second heating mechanism only applies when an addressing beam is used, and it results because the addressing 
beam necessarily Stark shifts the mobile sublevel. Therefore the trapping potential of an atom in that state is the 
sum of that due to the optical lattice light and the addressing light. The vibrational frequencies are thus different 
for the two hyperfine sublevels. Flip transitions between the storage and mobile sublevels can be made in two limits. 
If the pulse time is long compared to the inverse of the vibrational state splitting, so that the vibrational states are 
resolved, then the atom can make a transition to the new vibrational ground state, and there will be no heating. If 
the pulse time is shorter than the inverse of the vibrational state splitting, then the original vibrational wavefunction 
will be projected onto a superposition of states in the new basis. These levels will tend to dephase, and the atom 
will no longer be in the vibrational ground state when it is returned to the storage state. This heating is significantly 
reduced when the addressing beam is weakened, which requires that the microwave pulse be longer. Calculations of 
this heating effect have been performed in the context of making single qubit operations in a site addressable optical 
lattice, and the heating is ~ 10^^ vibrational energies per 30 fjs flip operation |27l |. 

Another heating mechanism is due to the photon recoil from the flipping radiation. In the case of microwaves and 
co-propagating Raman beams, the photon recoil is negligible. But for the orthogonal Raman beams, unless the pulse 
is slow enough to resolve the atom's vibrational states, the target atom will receive a significant photon recoil kick. 
In the Lamb-Dicke limit, the probability of vibrational excitation per pulse is not high, but this would still likely be 
the dominant source of heating in the compacting sequence. 

VI. SHIFT OPERATION 

The shift operator moves the lattice trapped atoms over a distance a/2, where a is the lattice constant. As mentioned 
earlier, all the atoms are initially prepared in a storage state, e.g. rriF = +1. It is easily seen from Eq. Q that the 
atoms will move when the linear polarization vector of one of the lattice beams is rotated. Rotating the relative angle 
of polarization ^ by tt moves all the atoms by a/2. Switching the internal state of the mobile atoms to the mobile 
state, and slowly returning the polarization to its original value, 9 = 0, returns the stationary atoms back to their 
original positions, while each mobile atom is moved forward another a/2, so that it is separated by a from its original 
location. Finally, the mobile atoms are returned to the storage state. 

As seen from Eq. (|2l, at = the atoms are confined only by the first term of the potential. As 6 increases, 
the second term starts to dominate the optical lattice potential. At 9 = tt/2, where the well-depth is Ui, only this 
state sensitive term confines the atoms. Further rotation to 9 = tt increases the well depth back to its original value, 
Uq. We assume that the relative angle between the polarization vectors of the lattice beams, 9, is a linear function 
of time, i.e. 9 = f3t, though it is conceivable that a non-linear variation may be useful to further reduce the heating 
under certain conditions. 

The highest vibrational heating occurs for the largest variation of the potential during the shift operation. For 
ground state Cs atoms with |mi?| = 1 in a 50 nm blue detuned optical lattice, the potential well-depth changes by 
a factor of approximately Uq/Ui = 8 during one shift operation. This ratio is higher in the case of the CO2 optical 
lattice, which is much further detuned from the atomic transition. The corresponding change in the frequency of the 
lowest vibrational levels is illustrated in Fig. O We now analyse the case of Uq = 100 iiK and a = 5.3 iim, in order 
to obtain an upper bound on the heating process. 

A. Vibrational heating 

Initially each lattice atom is in its vibrational ground state with energy Uq. The shift operation displaces the lattice 
potential and causes vibrational excitation of atoms. This heating can be quantified by the total energy of a moving 
particle relative to the initial state, < E > — Uq. We have employed both analytical and numerical techniques to get 
insight into the dynamics of atoms in a time- varying lattice potential. Application of the analytical approach is limited 
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FIG. 3: Change of the lowest vibrational eigenstates of the periodic potential for the atomic state \Cs 6S1/2, F = 4, mp ~ —1) 
as a function of the rotation of the relative polarization angle 6. The energy is measured relative to the potential minimum. 
The vibrational frequency for the transition — » 1 is 10.5 kHz for ^ = 0, and 3.66 kHz for — ^. 
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FIG. 4: The maximal vibrational heating of the initial ground state population AE as a function of the total timescale r of 
the shift operation. The analytical and numerical simulation results in the harmonic approximation are compared together 
and with the numerical simulation using the simplified periodic potential (see the text for explanation). The analytical and 
numerical results are indistinguishable on this scale. 



to the simplified case where both terms of the potential (01 are of the same magnitude. Numerical simulation was 
made using a Fourier grid representation of quantum states and operators and the Chebychev polynomial quantum 
propagation method. Details of the simulation methods can be found in Appendix 1X1 The explicit time-dependence 
of the potential was approximated in discrete steps, with a time step At I/Uq chosen to ensure the accuracy of the 
final results to within ^ 10^^. 



1. Uo = Ui case 



We first consider the simplified case where the two terms in Eq. (|2J) have the same magnitude, i.e. Uq = Ui. In the 
vicinity of the minimum, we can approximate the periodic potential of the lattice with a second order Taylor expansion, 
yielding the harmonic vibrational frequency uj — kyj2Uo/m, where k — 2it/\. Comparison with the results of direct 
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FIG. 5: The minimal fidelity F as a function of the shift operation timescale for the simplified case (see the text for explanation). 
The very good agreement between the analytical and numerical results proves that the classical picture of the process based on 
the coherent state representation correctly reproduces the atomic dynamics induced by the time-varying potential. For short 
timescales, the anharmonicity of the periodic potential slightly reduces the fidelity compared to the results obtained within the 
harmonic approximation. 

diagonalization of the periodic potential shows that the fractional difference from the true frequency is ~ 10^'^, and 
that the anharmonicity of the periodic potential (manifested in the deviation of uJk,k+i = iE^=k+i ~ Ey=k)/fi- from 
the harmonic frequency) is linear in the vibrational quantum number k and is only 3 % for k = 20. These facts 
justify the use of the harmonic approximation to get analytical insights into the process of vibrational heating. 

The motion of the simplified potential, induced by rotation of the polarization vector of one of the lattice beams, is 
linear in time and induces a transfer of the energy into the system E = rav^ = m^^. Here m is the atom mass, and 
Az is the total distance of the potential translation over the duration r. E is the most energy that can be transferred 
into vibrational motion during the beginning or end of a translation. Onset of the potential translation at time t = 
causes displacement of the initially stationary vibrational state with respect to the moving potential reference frame. 
The displacement transforms the stationary initial state into a vibrating coherent state with a maximal displacement 
of Xmax- This displacement is related to the transferred energy via E — ^mui'^x^^^, so x,nax — The maximal 

energy in the system at 6 — n is 2E. The analytical and numerical results for the harmonic and periodic potential 
show perfect agreement, as illustrated in Fig. ^ and hence justify the applied harmonic approximation and the 
coherent state representation. 

To further characterize the quality of the shift operation, we define a fidelity measure, F — = = 

Ol^"/)!^, which corresponds to the projection of the final coherent state onto the vibrational ground state of the 
translated potential. In this simplified case, the fidelity can be evaluated analytically, yielding F = exp{—^x'^^^) = 

ea;p(— As can be seen in Fig.[Sl the fidelity remains low for shift operation times less than 0.5 ms, whereafter 
it rises sharply. For shift times longer than 2 ms the fidelity is very close to unity, and thus provides a good estimate 
of the timescale necessary to preserve adiabaticity. 

2. Control of heating 

Exploration of this simplified model provides insight into the control and minimization of vibrational heating. 
The onset of the potential motion generates a coherent state from the initial vibrational ground state. In the moving 
reference frame, the potential motion first displaces the initial state to the negative momentum and negative coordinate 
region of the phase space, while setting its motion along the classical phase space trajectory corresponding to the 
coherent state. At the final instant of time when 9 = n, the actual motional phase of the coherent state accumulated 
during the evolution (e^*'^'^) determines whether stopping the potential will reduce or increase the total vibrational 
energy of a trapped atom. If the shift operation duration t is an integer multiple of the vibrational period, then the 
deceleration of the potential cancels the heating from the acceleration. It is thus possible to control and minimize the 
heating from this process. 
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FIG. 6; The maximal vibrational heating of the initial ground state population l\E, as a function of the total timescale r of 
the shift operation on the realistic potential (see the text for explanation). The heating has two contributions. The larger 
component originates from acceleration and deceleration of the potential during the motion. The smaller component comes 
from the phase of the vibrational motion of the coherent state at the point when the potential stops (9 = tt) and it is controllable 
by suitable timing of the shift operation as in the simplified analysis (see text). 



We next consider the realistic case where the depth of the potential depends on the relative linear polarizations of 
the lattice beams. For our parameters, the depth changes from C/q to L/q/S and back to Uq, when the angle between 
polarizations is rotated through tt radians. The vibrational frequency is therefore also modulated, accelerating the 
atoms when < tt/2 and deccelerating them when 9 > tt/2. This process of acceleration and deceleration of the 
motion, absent in the simplified analysis above, dramatically affects the final amount of vibrational heating of the 
atoms in the realistic potential. 

The real system can also heat due to the motional phase effect discussed in the simplified case. As in the simplified 
case, it can be eliminated here by suitable timing of the shift operation. Although for the case Ua/Ui ~ 8, it is a 
rather small portion of the total heating (see Fig. E)), it becomes significant when either the duration of the shift 
increases or the ratio Uq/Ui decreases (which ratio depends on lattice detuning and the choice of atomic states. 

The presence of two contributions to the heating energy is apparent from the fidelity plots shown in Fig[71 If the 
acceleration dominates the heating, the fidelity converges faster to its maximal value F = 1 compared to the case 
dominated by the controllable (oscillatory) contribution. 

Optimal conditions can be achieved by varying time duration r, well-depth Uq (see Fig. IS)) and the ratio Uo/Ui. 
When Uq/Ui = 8, the atoms can be shifted on a time scale of 5 — 6 ms without appreciable heating, as shown in 
Fig. O Moreover it can be seen that even for fairly fast shift operations, r < 2 ms, the atoms remain deeply trapped 
in the lattice, although recooling of the lattice may be required after many operations. The fidelity F provides a 
sensitive measure of error. Figure [T] shows that atoms usually stay in the ground state when r > 6 ms. 



Given that the initial lattice filling factor is only fi — 0.5, multiple shift operations will certainly be required to 
compact the lattice. Thus, it is important to ensure that the heating of the atoms is a well-behaved function of the 
number of shifts performed. In fact, it can be expected that the maximal heating scales as a linear function of the 
number of the shift operations, because of additivity of energy. There is no reason to expect a different scaling with 
the number of shifts if these are uncorrelated. When an enlargement of the computer size requires an increase of the 
number of the compacting operations beyond the tolerable heating limit, an intermediate recooling of the lattice by 
Raman sideband cooling can be applied to reset the ground state occupation to one. 



3. Realistic case 



4- Scaling 
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FIG. 7: Plots of the minimal fidelity F as a function of the shift operation timescale for the full potential with time varying 
ratio Uo/Ui. As this ratio decreases, the contribution to the vibrational heating from the acceleration of the potential motion 
is decreasing, and eventually becomes dominated by the oscillatory contribution. The latter has the same origin as that in the 
simplified analysis (no acceleration) and is fully controllable by suitable timing of the shift operation. The character of the 
fidelity plot in this case also shows slower convergence to F — 1 typical for the simplified case (Fig. |5J . This indicates that for 
durations r ~ 3 ms, a higher fidelity can be obtained with the potential Uo/U\ = 3 than for U^jUx — 2. 
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FIG. 8: The minimal fidelity _F as a function of the well depth [/q for the realistic case with Uo/U\ = 8. 



VII. EFFICIENCY OF COMPACTING 



To systematically analyze the efficiency of compacting, we in turn consider finite 1-D, 2-D, or 3-D lattices with 
50% initial site occupation Pocc- From imaging during laser cooling, we have a map of the lattice occupation. Our 
goal is to move the atoms so that they form a single contiguous block with no vacancies, i.e., a perfect finite lattice. 
The primitive operation in this analysis is the shift of an atom or group of atoms by one lattice site. We define the 
cost (time) for this operation to be one unit. The physical implementation of this operation has been discussed in 
section IVn where we showed that moving a group of n atoms through one lattice site requires two movements of the 
trapping potential (2 elementary shifts) and 2n flip operations (see Fig. Note that there are only two elementary 
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shifts because all atoms move in a shift. For fauh tolerant computation |23| we need to have 0{n) parallel operations 
|29] , so for evaluating the scalability of the scheme we assume that the spin flips can be achieved in parallel. This is 
a good assumption, since the time required for a spin flip is much smaller (by a factor of 100-1000) than the time to 
move the trapping potential. In fact, we can ignore the spin flip cost altogether. We also ignore the cost of classical 
computations required to plan and implement the compacting procedure since, with proper optimization, this cost 
scales at most linearly with the total number of lattice sites. Therefore, the total cost is essentially determined by 
the shift operations. 

We first study the 1-D lattice to show the solution in a simple setting. Then we consider the 2-D lattice, which 
presents additional challenges. Finally, we generalize the techniques of the 2-D lattice to the 3-D lattice. 

A. One Dimensional Lattice 

The one dimensional lattice i„ has n sites. The site occupation probability is Pocc- Our goal is to move the atoms 
so that they form a line of atoms with no gaps in between. 

Here, we use the obvious algorithm, which we call COMPACT, to remove the vacancies. The algorithm moves 
the atoms to the left so that all the vacancies move to the right. Thus after running the algorithm we end up with a 
line of atoms at the left of the lattice. 
COMPACT: 

1. Find the leftmost vacancy v. If there are no atoms to the right of this vacancy, STOP. 

2. Let G be the set of atoms that are to the right of v. Shift all the atoms in G left by one step. GOTO^ 

We say that a vacancy is at the right side if all the atoms are to its left. Since there are at most n vacancies and 
each operation can take one vacancy to the right side, the number of operations needed is n. On average there are 
(1 — Pocc)n vacancies, so the expected number of shifts required is (1 — Pocc)n. Also, if there are exactly v vacancies 
that are not on the right side, then the number of operations needed will be v. 

Clearly, the cost of this algorithm in the average case can be lowered by finding the center of mass of the set of 
atoms and compacting them around this. However, this procedure does not improve the worst case cost. 

B. Two Dimensional Lattice 

We consider a finite two dimensional square lattice L„.„ with N ~ sites. We denote a sublattice as Sij where 
(1 < i < J < n). Si J consists of all rows between and including rows i and j. Let N{Sij) denote the total number of 
atoms in the sublattice Sij. Let represent the number of atoms in the row k. 

In this case the compacting procedure has two parts: 

1. ^BALANCE (5i.„): The atoms are first moved so that they are equally distributed among all the rows. The 
general idea is to divide the lattice into two halves, e.g. the top half and the bottom half (Fig. (SJ. Then, atoms 
are moved from the top half to the bottom half or vice versa, so that each half contains the same number of 
atoms (a difference of one atom is allowed when the total number of atoms is odd). After this, each half of the 
lattice is given as input to this procedure again, thus recursively balancing all the rows. Fig. |3 illustrates this 
procedure for a 7x7 lattice. A more detailed description of the algorithm can be found in appendix ^ The 
recursive procedure can be summarized graphically by a binary tree structure (Fig. llOfl . as described below, in 
order to find its cost. 

2. ^ROW-COMPACT: The atoms in each row are then compacted as in the 1-D case, except here all rows can be 
compacted in parallel since mobile atoms can be moved together as a group. Fig. illustrates this procedure 
for a 7x7 lattice that has already been balanced according to Fig. |21 

Now we analyze the cost of the entire procedure. First we analyze the cost of ROW-COMPACT. It is the same as 
doing the 1-D compacting on each row, but we can do all the rows in parallel so the worst case cost is n jSOl (there 
being only n sites per row, thus requiring at most n — 1 moves to get an atom to the left). The average case cost is 
again (1 — Pocc) since there are on average (1 —poccin vacancies per row. Additionally, if v is the maximum number 
of internal vacancies amongst all rows, then the cost is v since we have to move all the vacancies beyond the last atom 
on the right. 

Let us now consider the cost of the balancing procedure. We can represent the recursion by a binary tree of depth 
[log2 n] as shown in Fig. ^| The nodes of the tree are labelled by the first and last rows of the sublattice Si,j passed 
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FIG. 9; Here we show the movement of atoms during the balancing process for a 7x7 lattice. Figure (a) shows the first halving 
and balancing, corresponding to the 0th level of the binary tree from Fig. 1101 Figure (b) shows the balancing corresponding 
to the second level of the tree. Figure (c) shows the balancing corresponding to the third level of the tree. Figure (d) is the 
final result where all the rows are balanced. The horizontal shaded rectangles show the partitioning of the lattice into halves. 
The grey circles represent atoms that are to be moved in the corresponding step. The black circles represent atoms that do 
not move in that step. 



Depth/Level 




FIG. 10: The figure shows the binary tree describing the recursion of BALANCE when it is given Si,7 i.e. a 7x7 lattice viz. 
the one considered in Fig. |U] The nodes of the tree show the starting and ending rows of the sublattices balanced in each step. 
The root of the tree has depth or level 0. The maximum depth is [logj n] which is 3 in this case, since n = 7. Note that at 
the last level, i.e. at the leaves of the tree, all the lattices consist of one row only. Hence no balancing needs to be done at this 
point. Thus the total number of balancing steps is [log2 n] . 



as input to BALANCE. The set of nodes at the same depth d are said to be at level d. The root is defined to be at 
depth 0. During each recursive step, the number of rows in the sublattice is halved, so the depth of the tree is [log2 n] 
[31|. The total number of shifts required is the sum of the number of shifts required at each level of the tree. Notice 
that the number of atoms to be moved during a balancing step is not a problem since they can move in parallel as a 
group. We basically need to count the number of shifts that the atoms have to make at each level. For a sublattice of 
k rows the atoms only need to make \k/2] shifts at most, since that is the largest number of rows of the two halves 
that are being balanced. At depth d the number of rows in a sublattice is at most 23tt 113 ■ There are 2''+-'^ such 
sublattices to be balanced at depth d. Now, the balancing of all these sublattices can be done in parallel: first move 
all the atoms in all these sublattices that need to be moved up. This takes at most shifts. Then shift all the 
atoms in all the sublattices that need to be moved down. This also takes at most ^^tt shifts. Thus the total number 
of shifts required is ^ + + j...{\\og2 n]tcrms)) < f + § X]^^Lo(5)" ~ where the first term in the summation is 
only n/2 (not doubled) because there is only one direction to move the atoms the first time the lattice is halved. Also 
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FIG. 11: Here we show the movement of atoms during the row-compacting process for the 7x7 balanced lattice obtained at the 
end of the balancing shown in Fig. |5| The grey circles represent atoms that are to be moved in the corresponding step. The 
black circles represent atoms that do not move in that step. 

note that one would expect to be able to balance each level of the binary tree in one step on average. This would 
imply that the average complexity of balancing is « logj n shift operations. Our preliminary simulations suggest that 
this is indeed the case. 

Putting together the above component of the cost analysis we see that in the two dimensional lattice with N = 
sites, the compacting procedure takes at most n+ = ~ f steps (ignoring terms 0(log A^)). 

C. Three Dimensional Lattice 

We consider a three dimensional cubic lattice Ln,n.n with N = rr' sites. The compacting procedure for this three 
dimensional case is a generalization of that for the two dimensional case. One key difference is that the notion of rows 
needs to be defined. Say, the lattice is in the positive octant of the coordinate system, and its three sides coincide with 
the axes. Taking the lattice spacing to be the unit along the axes, the lattice sites are given by the coordinates (i, j, k) 
with 1 < i, j, k < n, where the first coordinate is the x coordinate, the second coordinate is the y coordinate and the 
third coordinate is the z coordinate. We define row Rij to be the ordered set of lattice sites | 1 < / < n}. 

The Z*'' position in the row Rij corresponds to the site 

Define a sublattice Sij-k,i with {l<i<j<n,l<k<l<n)to consist of the set of lattice sites {(a, 6, c) | 1 < 
a<n, i<b<j;k<c<l}. Let N{Sij:k.i) denote the total number of atoms in the sublattice Sij-k,i- 

The lattice compacting procedure consists of 

1. ■^BALANCE(S'i^„;i^„,0): The atoms are moved so that they are equally distributed amongst the rows. The 
second parameter in the function call corresponds to the current recursion depth. The general idea is the same as 
in the case of the two dimensional lattice. The difference is now that we alternately balance the two halves made 
by planes parallel to the xy plane and the xz plane, using the recursion depth d to implement this alternation 
of balancing. 

2. ■^ROW-COMPACT: The atoms in each row are then compacted in a fashion similar to the 1-D case, except now 
the operations for all the rows can be done in parallel since we can move a three dimensional group of atoms 
together in a single shift. 

Now we analyze the cost of the entire procedure. First we analyze ROW-COMPACT. The analysis is identical to 
that for the two dimensional case. For the sake of completeness we mention the results here. The worst case cost is 
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(xy-balance) 

1 (xz-balance) 

2 (xy-balance) 

3 (xz-balance) 

4 (stop) 




FIG. 12: The binary tree describing the recursion of BALANCE when it is given S'i,3;i.3 as input, i.e. n = 3. The nodes of 
the tree show the indices of the sublattices balanced in each step. The root of the tree has depth or level 0. The maximum 
depth is 2[log2 n] which is 4 in this case since n = 3 to the even levels correspond to balancing w.r.t. to the first set of indices 
(halving plane is parallel to xy plane), while the odd levels correspond to balancing w.r.t. the second set of indices (halving 
plane is parallel to the xz plane). Note that at the last level, i.e. at the leaves of the tree, all the lattices consist of one row 
only hence no balancing needs to be done. Thus the number of non-trivial balancing steps is 2 [log2 n] . 

n. The average case cost is (1 — Poccjn-. Additionally, if v is the maximum number of internal vacancies amongst all 
rows then the cost is v. 

The analysis of balancing for the three dimensional lattice is similar to the two dimensional case. During each 
recursion step the lattice is halved along the xy or xz plane so the depth of the recursion is 2 [log2 n] . We can 
represent the recursion by a binary tree of depth 2 [logj n \ , where even levels correspond to halving the sublattice 
parallel to the xy plane and odd levels correspond to halving the sublattice parallel to the xz plane. The sum of the 
number of shifts required for each level of the tree gives us the total number of shifts required during the execution 
of the algorithm. We will do the counting separately for the odd and even levels. Since mobile atoms can be moved 
in parallel, we only need to count the number of shifts that the atoms have to make during each balancing step. 
For a sublattice of k rows along the direction in which we halve the sublattice atoms, we only need to make \k/2~\ 
shifts. At depth 2d and 2d + 1 the relevant number of rows in the sublattice is 23tt (taking the root to be at depth 
0) [S^. As in the 2-D case, we parallelize the required shifts for all the sublattices at a given depth k, by doing 
all the shifts in the same direction in parallel. Since there are two directions, the number of shifts required at each 
depth is twice that required for a single sublattice at that depth. Then the total number of shifts at the even levels 
is n/2 + n/2 + n/4-|-n/8...([log2nlterms)) < n/2 + n/2Y,n=oi'^/^y'' = 3/2n. The first term in the sum is only n/2 
(not doubled) because there is only one direction to move the atoms the first time the lattice is halved. For the odd 
levels a similar argument holds, except that the first term is n because there are two sublattices at that stage and 
balancing them may require moving atoms in opposite directions. Then the total number of moves required at the 
odd levels is n + n/2 + n/4 + n/8...([log2 n]terms)) < «X]^o(-'^/^)" ~ Adding together the total moves at the 
even and odd levels gives us the total number of moves required: 3/2n + 2n — 7/2n. Thus, in the three dimensional 
lattice with N — sites, the lattice can be compacted in at most 7/2n + n — 9/2n — %/2\/N steps (ignoring terms 
0(log7V)). As in the case of the two dimensional lattice, one would expect that the balancing at each level of the 
binary tree can be typically done in one step. This would imply that the average cost of balancing is 2[log2]?T. steps. 
However, we do not have a rigorous proof of this estimate. 

VIII. CONCLUSION 

In conclusion, we have presented an experimentally viable scheme to perfectly initialize a quantum computer based 
on neutral atoms trapped in an optical lattice of large lattice constant. The proposed compacting scheme allows 
preparation of a perfect optical lattice of orthorhombic structure, where each lattice site is occupied with a single 
atom in a specific internal level and motional ground state. We have proposed physical mechanisms for realization 
of the two elementary compacting operations, the flip of the internal state of atoms and their shift within the lattice 
structure, and have analyzed their efficiency in some detail. Particular attention was devoted to the study of the 
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motional heating during the flip and shift operation. Mechanisms to control this heating were proposed that are 
based on the analytical and numerical solutions of the atomic dynamics on the time-dependent lattice potential. 
Our analysis of the complexity of the compacting process demonstrates its scalability with the size of the quantum 
computer. 

We point out that all of the components discussed here are feasible with current technology. Our results show that 
we can achieve a perfect optical lattice in less than one second starting from a half filled lattice of approximately 8000 
sites. The procedure can also be used to complement other proposals that provide near complete but not perfect filling, 
provided they are carried out in a site addressable optical lattice. Furthermore, the lattice compacting techniques 
proposed here can be adapted to the execution of quantum logic gates in site-addressable optical lattices. The control 
of individual atoms in an optical lattice envisioned here would set the stage for the construction of a neutral atom 
quantum computer. 
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APPENDIX A: METHOD OF NUMERICAL SIMULATION 

In this appendix, we summarize the two parts of our approach to numerical simulation of atomic motion on 
an optical potential. First, we introduce the Fourier grid method, a discretization of the Schrodinger equation 
suitable for computational implementation. Second, an accurate and efficient approach, based on a Chebychev 
polynomial expansion of the time evolution operator, is presented. In both paragraphs, we discuss the advantages 
and disadvantages of the methods. 

1. Fast Fourier grid method 

A combination of continuous and discrete degrees of freedom characterize full scale dynamics of most quantum 
mechanical system. The discrete degrees, at least in cases of smaller dimensionality, are often easily handled, but the 
continous degrees of freedom require that some form of discretization scheme be used. 

Finite differences introduced grid methods into quantum molecular dynamics simulations, but, since they suffered 
from poor accuracy and slow convergence, they were soon replaced by other methods. The Fourier method does 
not suffer from the same shortcomings, and allows for rapid convergence and efRcient sampling of system phase 
space. 34, 35] 

The Fourier basis functions, 9k = e^'^^^^l^ ^ form an orthogonal basis, and can be used to represent a system 
wavefunction or operator to arbitrary accuracy. 

Af-l 

\^{x)) = GkOk (Al) 

k=-N 



2N 

1=1 

For any non-local system operator, e.g. the kinetic energy operator in the position representation, the evaluation 
of the operator, O, scales as the square of the number of grid points. In contrast, the potential energy operator, 
U{x), which is local in position space, scales linearly with the number of grid points and is thus simple to evaluate. 
The advantage of the Fourier method is that the kinetic energy operator, T{p), can also be evaluated as a diagonal 
operator if the wavefunction is expressed in the momentum basis. Consequently, the action of a given Hamiltonian, 
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H{x,p) = T{p) + U{x), can be evaluated as 

H^{x)^— e'-^T{p) J2 e~"-^^{x) + U{x)^{x) (A3) 

The first term describes the action of the kinetic energy operator. It expHcitly takes into account the two transforma- 
tion steps: i) the Fourier transform into the momentum representation followed by application of the diagonal kinetic 
energy operator; and ii) the inverse transform back into the position representation. Additonally, to facilitate a direct 
application to physical problems, the Fourier transform is defined on the finite interval [—L/2, L/2], giving rise to the 
2TT / L term in the exponent. 

The Fourier method is limited to bandwidth limited functions that decay exponentially fast as the argument 
approaches infinity. For such functions the Whittaker-Kotel'nikov-Shannon theorem ensures that, with a finite number 
of sampling points, the value of the function can be determined to any desired accuracy. 36] In the case of time- 
dependent problems, the choice of sampling points has to ensure that the wavefunction remains confined at all time 
to the space defined by these points. 



2. Chebychev propagation method 

The Chebychev propagation algorithm employs the orthogonal Chebychev polynomials to approximate the quan- 
tum mechanical time evolution operator, U{t), or any other operator functional. In contrast to other interpolation 
methods, the Chebychev method does not suffer from inherent instabilities that limit the order of the expansion, and 
consequently also the attainable accuracy. Moreover, the interpolation error is always uniform, and can, due to 
the possibility of very high order expansions, be reduced to the order of the computer's numerical precision. 

Because Chebychev polynomials are defined on the interval [—1, 1], it is necessary to map the Hamiltonian onto 
this interval using the linear transformation 

H' - 2 ^'-^"'"-^ _ /. (A4) 

Here H is the Hamiltonian of the system, Emax and Emm are the estimates of the largest and smallest energy 
eigenvalues of H, and / is the identity operator. For the next step of the Chebychev propagation algorithm, the 
Chebychev expansion coeffecients are calculated: 




Here bk is the expansion coefficient for the Chebychev polynomial of order k, Tk{H'). d{At) is a global phase factor, 
6{At) = exp{—{i/h){AE/2 + Emin)At), arising from the mapping step. The appearance of the Bessel function, Jfe, is 
attractive since it provides an automatic cut off point for the expansion. As the order of the expansion increases beyond 
the size of the argument, , the associated expansion coeffecients will approach zero exponentially fast, The 

time evolution operator can now be approximated as 

N 

u{Aty = ^^^^ (^') (A^) 

fe=0 

Here hk is the expansion coefficient for the k'th order Chebychev polynomial, Tk{H'). More efficient than calculating 
and applying the full expansion to the system wavefunction, however, is recursively updating the system state. This 
recursive approach is particularly attractive since it permits the expansion to be terminated when the desired precision 
has been attained. Recursion is carried out as follows: 



\^{t + At)), = \^{t)) 
\^{t + At))^ = H'\^{t)) 

\^{t + At))k - 2H'\^{t + At))k-i-\^{t + At))k-2- 



(A8) 
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APPENDIX B: 2-D BALANCE 

2BALANCE(5^,™) see Sec. |VIl|for definitions. 
INPUT: S^j 

GOAL: Balance the rows in Si_j so that V k, I we have \nk — ni\ < 1, where i < k,l < j . Equivalently, each row has 
at least rimin — — * + atoms. 

ALGORITHM 

1. Let Z = (j - i + 1) denote the number of rows in Si^j. Ul = 1 RETURN. 

2. Let m = i + [l/2\ be a middle row of Sij. For balancing, we need to have rireq — (jn — i + l)nm,in atoms in the 
sublattice Si^m- 

3. If {N{Si^rn) > rireq), shift {N{Si^rn) - flreq) atomS down frOm Si^rn tO Sm+1,] 

4. If iN{Si^,n) < rireq), shift {Ureq - N{Si^„i)) atOmS up from S„i+l,j to Si^jn- 

5. 2bALANCE(S'„,+ij). 

APPENDIX C: 3D-BALANCE 

3BALANCE(5j,™;fc,i,d+ 1) see Sec. EHlfor definitions. 
INPUT: Sublattice Si,j-k,i and recursion depth d. 

GOAL: Balance the rows in Sij-k,i so that Va.6 such that i < a,b < j and Vc.d such that k < c,d< I we have 
\na,c — "-b,d| < 1, where Ua.c represents the number of atoms in the row i?a,c- Equivalently, each row should have at 

least Ujrtin = L (j-jyii^Ci-fc+i) J atoms. 
ALGORITHM: 

1. If i = j and fc = / then RETURN. 
Only one row so no balancing needed. 

2. If d is odd, gotoHnielse if i = j 3BALANCE(S'jj.fcj, + 1). RETURN. 

3. If sublattice Sij-k,i contains only one xy plane (i.e. i = j ), then gotolTUI 

4. Let / = J — i + 1 be the number of planes parallel to the xy plane. Then m = i + \J-/'2\ defines a middle plane 
parallel to xy plane. 

5. Let rireq = nminiTn — i + 1)(Z — + 1) be the minimum number of atoms that are required to be in Si^rn:k,i- 

6. If {N{S.i^rn-k,i) > rireq), shift (N(Si^m-k,i) - rireq) atoms fr 

7. If (iV(S'i,m;fe,i) < rireq), shift Ureq " N{St^rn;k,l) atOmS frOm Sm+l,j;k,l tO Si^rn;k,l- 

8. 3BALANCE(5™+i,,;fe,z,d+l) . 

9. RETURN 

10. If fc = ; ^BALANCE{Sij,k,i,d+ 1). RETURN. 

11. Do steps0|to|Hl replacing indices with indices {k,l). 

12. RETURN 
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